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ABSTRACT 

We present the results of a search through the photometric database of eclipsing Kepler binaries (Prsa et 
al. 2011; Slawson et al. 2011) looking for evidence of hierarchical triple star systems. The presence of a 
third star orbiting the binary can be inferred based on binary eclipse timing variations. We apply a simple 
algorithm in an automated determination of the eclipse times for all 2157 binaries. The "calculated" eclipse 
times, based on a constant period (Keplerian, or two-body) model, are subtracted from those observed. The 
resulting O-C (observed minus calculated times) curves are then visually inspected for periodic residuals in 
order to find triple-star candidates. After eliminating false positives due to the beat frequency between the 
~ 1/2-hour Kepler cadence and the binary period, 39 candidate triple systems were found. The periodic O-C 
curves for these candidates were then fit for contributions from both the classical Roemer delay and so-called 
"physical" delay, in an attempt to extract a number of the system parameters of the triple. We discuss the 
limitations of the information that can be inferred from these O -C curves without further supplemental input, 
e.g., ground-based spectroscopy. Based on the limited range of periods for the triples to which this search is 
sensitive, we can extrapolate to estimate that at least 20% of all close binaries have tertiary companions. 
Subject headings: stars: binaries: general — stars: formation — stars: triple — stars: 



1. INTRODUCTION 

Triple star systems are appealing objects for study for a 
number of reasons. The orbital architecture and masses of 
the constituent stars can inform us about the not-so-well un- 
derstood process of the formation of systems of multiple stars 
(see, e.g., Boss 1991; 1995; Bodenheimer et al. 2000; Sterzik, 
Tokovinin, & Shatsky 2003; Bate 2009; Reipurth & Mikkola 
2012). As one example, it is known that close binary sys- 
tems cannot have formed in their current configurations; dur- 
ing their protostellar phase the stellar radii would have been 
much too large to fit inside their current orbits. The presence 
of an orbiting third star in the system could provide a natu- 
ral mechanism, through Kozai cycles (Kozai 1962) with tidal 
friction, for the initially wide binary to lose angular momen- 
tum and become close (Kiseleva, Eggleton, & Mikkola 1998; 
Eggleton & Kiseleva-Eggleton 2001; Fabry cky & Tremaine 
2007). This mechanism has also been proposed as a way to 
explain the blue-straggler stars found predominantly in glob- 
ular clusters (Perets & Fabry cky 2009). The orbital architec- 
ture of a triple system can also in principle inform us about the 
final contraction of the interstellar cloud that formed the sys- 
tem, provided the dynamical evolution of the system has left 
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the initial configuration relatively unaltered (see, e.g., Boss 
1991; Bate 2008; Reipurth & Mikkola 2012). 

Moreover, understanding the relative frequency of binaries 
vs. triples and quadruples (see, e.g., Tokovinin et al. 2006; 
Pribulla & Rucinski 2006; Raghavan et al. 2010) is important 
in anticipating what other unseen stars in any particular sys- 
tem may be present. The hypothetical presence of such bod- 
ies may be important in explaining various effects that are ob- 
served in these binaries, but not otherwise explained (see, e.g., 
Eggleton & Kiseleva-Eggleton 2001, and references therein). 
Finally, while studies of binary star evolution, and especially 
the phases involving mass transfer, have dramatically trans- 
formed our overall understanding of stellar evolution and the 
exotic remnants, such as binary neutron stars, that are left in 
the late phases, studies of the little-explored triple star evolu- 
tion promise to involve yet several more layers of complexity. 

There are at least five ways of finding triple-star systems. 
These include (i) visually resolving bound star systems, in- 
cluding with adaptive optics and optical/IR interferometry 
(see, e.g., Tokovinin et al. 2006; Rucinski, Pribulla, & van 
Kerkwijk 2007; Raghavan et al. 2010). (ii) Observing the 
presence of three different stellar spectra in an apparently 
single object provides an excellent starting point for the dis- 
covery of triples (see, e.g., Zucker, Torres, & Mazeh 1995; 
D'Angelo, van Kerkwijk, & Rucinski 2006). (iii) Doppler 
spectroscopy carried out over periods at least as long as the 
binary period in the system, and a substantial portion of the 
triple period, is the most informative (see, e.g., Carter et 
al. 201 1). (iv) Direct observations of eclipses by all three bod- 
ies is also exceptionally interesting, but such systems are rel- 
atively rare (see, e.g., Carter et al. 201 1; Derekas et al. 201 1; 
Carter et al. 2013). Finally, as has been done for more than 
a century (v) long-term timing of binary eclipses can reveal 
periodic perturbations to the otherwise linear progression of 
eclipse times with cycle number (see, e.g., Irwin 1952; Fab- 
rycky 2010; Steffen et al. 201 1; Gies et al. 2012; Borkovits et 
al. 2013). It is the latter approach which is the subject of this 
paper. We also note that this method of timing variations has 
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been used to great success in measuring orbits and masses of 
multi-planet systems (see, e.g., Holman et al. 2010; Lissauer 
et al. 2011; Carter et al. 2012), though the mass and period 
ratios of the perturbers are different in planetary systems vs. 
triple star systems. 

However, each of these methods suffers from some limi- 
tations, and each probes different regimes in the ratio of the 
binary period to that of the triple. In the case of timing binary 
eclipses, this can be done quite accurately from ground-based 
measurements, at least on bright objects, and such studies 
have provided substantial hints of the presence of third bodies 
(see, e.g., Pribulla & Rucinski 2006). The difficulty here has 
been that ground-based eclipse timing studies are subject to 
frequent interruptions due to the diurnal, lunar, and seasonal 
cycles, not to mention the weather. In this work we make use 
of three years of nearly continuous observations by Kepler of 
some 2000 eclipsing binaries to identify candidates for triple 
star systems. 

The Kepler mission (Borucki et al. 2010; Koch et al. 2010; 
Caldwell et al. 2010) has been observing some 157,000 stars, 
including ^2000 eclipsing binaries, for the past three years. 
The continuous monitoring of these eclipsing systems, in 
combination with the exquisite high photometric precision 
of the Kepler mission (Jenkins et al. 2010a; 2010b), is un- 
precedented in the history of observational astronomy. As a 
result, this photometric data set of eclipsing binaries is able 
to make a serious contribution to the endeavor of identifying 
promising triple candidates for followup studies with Doppler 
spectroscopy. Already, the Kepler observations have yielded 
some five triple star systems identified directly by third-body 
eclipses of the binary (Carter et al. 201 1; Derekas et al. 201 1; 
Slawson et al. 2011) while a number of others have been in- 
ferred to be triples by evidence for systematic eclipse timing 
variations ("ETVs") of binaries (Fabrycky 2010; Slawson et 
al. 201 1; Steffen et al. 201 1; Carter et al. 2013). The Slawson 
et al. (2011) catalog of binaries, in which ten of these triples 
are briefly mentioned, was based on only 120 days of Kepler 
data, whereas approximately an order of magnitude more data 
now exist. 

In this study we present the results of a comprehensive 
search of the Kepler data base of binary systems for evidence 
of the presence of a third star. This was done by searching 
for periodic features in so-called O-C curves (observed mi- 
nus calculated eclipse times) of some 2000 eclipsing binaries. 
We find 39 good candidates for triple stars. In addition to ex- 
hibiting the periodic variations in the O-C curves indicative 
of a triple system, several of our candidates feature additional 
evidence for being triple. For example, two of the systems 
have third-body eclipses, while seven of them exhibit secu- 
lar variations in the depths of the binary eclipses indicative of 
a precessing orbital plane of the binary. As we show, 19 of 
the systems exhibit dominant classical Roemer delays, while 
another 1 1 have dominant physical delays (due to perturba- 
tions to the binary "clock", i.e., its orbital eclipse period). 
The especially interesting feature of these candidates is that 
we can directly follow perturbations to the binary orbit and/or 
the classical Roemer delay continuously over several cycles 
of the triple. 

The processing of the Kepler data for the 2157 eclipsing 
binaries is described in |2] Production of an O-C curve 
for each system is discussed in ^3] while an overview of our 
triple-star candidates is presentedin § |3.3| Expressions for the 
various effects that appear in the O-C curves are given quan- 
titatively in Our approach to the analysis of the O-C 



curves, in order to extract as much information about the 
physical system parameters as possible, is described in #5] 
Our results for the 39 triples found in the search are presented 
in ^6] We discuss the limitations on the determination of 
system parameters using only the Kepler eclipse timing data, 
without supplemental information that could be provided by 
ground-based spectral observations (and in some cases by the 
Kepler data themselves). All of these systems will require 
such follow-up observations in order to definitively determine 
the masses of the three stars and the orbital elements. In 
we discuss our results, with emphasis on what can be learned 
from only the O-C curves. Finally, we attempt to estimate the 
fraction of close binaries with tertiary stars of orbital periods 
< few years. 

2. DATA PREPARATION 

2.1. Kepler binary data set 

The data we use for this study are long-cadence (LC) 
lightcurves for all binaries published in the latest Kepler 
eclipsing binary catalog (Slawson et al. 2011). We used all 
the files from Quarter 1 through Quarter 1 3 which were avail- 
able for retrieval from the Multimission Archive at STScI 
(MAST). The data used had all been reprocessed with the 
PDC-MAP algorithm (Stumpe et al. 2012; Smith et al. 2012), 
which removes much of the instrumental noise from the flux 
time series while retaining the bulk of the astrophysical vari- 
ability in sources. For each quarter, we normalized the flux 
series to its median value, and then stitched the quarters to- 
gether into a single file for each source. 

2.2. Filtering the data 

The next step in the data processing was to apply a high- 
pass filter, based on the known period of the binary system. 
We took the stitched 13 quarters of data, described in section 
2.1, and filtered out the low frequencies (starspot activity, in 
particular), in the following way. First, the data were con- 
volved with a boxcar function of duration equal to the known 
binary period. Second, the smoothed data were subtracted 
from the unsmoothed data. Frequency components below the 
frequency of the binary orbit are thereby largely removed, 
while leaving temporal structures that are shorter than the bi- 
nary orbital period. The eclipses themselves are essentially 
unaffected. 

The reference epoch for all times in this paper is Barycen- 
tric Julian Day 2454900. 

3 . ECLIPSE TIMING ANALYSIS :0-C CURVES 

3.1. Measuring Eclipse Times 

The baseline algorithm we utilized for determining the 
eclipse times consists simply of testing each flux point in the 
Kepler data set for a local minimum and fitting a parabola to 
the lowest three points in the local minimum. Then the fit- 
ted parabola is used to interpolate between Kepler samples to 
find a more accurate time of eclipse minimum. As we show, 
this algorithm is quite good for short orbital period binaries, 
but begins to lose accuracy for longer-period binaries when 
the eclipse duration may consist of a substantial number of 
Kepler long-cadence samples. To carry out our initial search 
for periodic variations in the O-C curves, we used this ba- 
sic algorithm exclusively. However, after interesting systems 
were identified, we recalculated more accurate O-C curves 
using a better algorithm that involves more of the eclipse pro- 
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file (T. Borkovits, unpublished) for a handful of the binaries 
with periods with Pbin ^ 6 dayaj 
The parabola to be fit is of the form: 

F„ = a(t„- Stf +F min (1) 

where n = 1,2, or 3; ?i = -1, = 0, and ?3 = +1; and St is 
the offset of the time of the minimum with respect to the time 
of the point with the lowest flux of the three Kepler samples. 
The times are all dimensionless, and are in units of A?lc = 
1765.46 sec, the Kepler long-cadence sampling interval. We 
note that the parameter a in this expression implicitly encom- 
passes information about the relative sizes of the stars, limb 
darkening, orbital inclination, and so forth. Presumably for a 
given binary system this parameter remains a constant, though 
in practice, effects such as time-varying starspots, can slightly 
modify a. 

Since not all binary eclipses are well represented by a sim- 
ple quadratic function near minimum, we also considered a 
quartic shape. This is the next simplest shape for any sym- 
metric eclipse profile. Because there are four parameters that 
describe a symmetric quartic, this would require four or more 
flux points to fit. Five is the minimum number of points 
in a symmetric arrangement which can have a lowest flux 
point with two higher-flux points on either side. However, we 
judged this to be too many to use for the shortest period bina- 
ries - in some cases, the eclipse is only a few Kepler cadence 
points wide. Thus, to get a flavor for how a quartic might fit, 
we utilized a function of the following form: 

F„ = a(t n - Stf + pa(t n - St) 4 +F min (2) 

where the parameter /? was fixed at a representative value of 
0.3. Thus, there are still only three parameters to fit analyti- 
cally to three data points. Again, note that all the times are di- 
mensionless (i.e., in units of A?lc)- We also tried other values 
for f3, but found no improvement (i.e., reduced rms scatter) in 
the "quartic" algorithm. 

Once we found a potential eclipse time, and a correspond- 
ing value of F m i n we required that it be less than a certain 
threshold flux in order to be judged an actual eclipse and not 
just an uninteresting local minimum in the flux. Formally, we 
somewhat arbitrarily required that 

F mta < 0.4 -Fed +0.6 (3) 

where F ec i is the flux at the bottom of the primary eclipse in the 
folded light curve, and recall that the fluxes are all normalized 
to unity. In some cases, this allowed the secondary eclipse to 
also be picked up, but these were distinguished by the ~T80° 
phase shift from the primary eclipse. 

In general, the quadratic function produced better results 
than the quartic, i.e., less scatter in the O-C curves, but 
yielded a comparable number of candidate triple stars. Both 
funct ions were equally susceptible to spurious periodicities 
(see g3T2]>. 

8 After this work had essentially been completed, we developed a more 
sophisticated eclipse timing code that formally cross-correlates the epoch- 
folded binary light curve with the Kepler data train. We found all of 39 of the 
triple-star candidates with this improved code, including four new candidates 
that the original search missed. The quality of the O — C curves was hardly 
changed for most the systems with I\ m < 10 days, but there were some im- 
provements, i.e., lower scatter, for a few of the longer period systems. In eight 
cases, where the O—C curve significantly improved over the simple quadratic 
fitting algorithm, and where the O-C curve had not already been upgraded 
using the Borkovits (unpublished) code, we used those O — C results rather 
than the original. 



As a separate piece of the analysis, we also deliberately 
found the times of the secondary eclipses. However, in this 
work we do not directly utilize their O-C curves in the tim- 
ing analyses. We do discuss what supplemental information 
the secondary eclipses can yield in the case of eccentric bina- 
ries. We also tabulate which systems have secondary eclipses 
whose O-C curves exhibit different behavior than that of the 
primary eclipse. 

Finally, we note that even though the nominal separation of 
the flux points in the long-cadence mode, A?lc, is 1765.46 
sec, we were able to determine the times of eclipse minima to 
a typical empirically determined accuracy of ^20-100 sec, 
or < 5% of the timing metric. We list the rms residuals to the 
model fits for each source among our tabulated results. 

3.2. Searching for Interesting O — C Curves 

As we search for potential triple-star signatures among the 
O-C curves, we find many that exhibit spurious periodicities. 
These false positives are most often due to a beat between 
the frequency of the Kepler cadence and the frequency of the 
binary orbit. The two prominent beat frequencies are given 
by: 

/beat, i = fuc ~ /bin • int ( 4^ J ( 4 ) 

V/bin/ 

V/bin/ 

where "int" gives the truncated integer value, and /},;„ = 1 /fW 
and /lc = 1 /AfLc- For each O-C curve that we compute, we 
display these two prominent expected beat periods. If there 
is a match between a predicted beat period and the detected 
period in the O-C curve, that object is eliminated as a possi- 
ble triple-star candidate. We note that these beat frequencies 
change (sometimes fairly obviously) during the course of a 
year. This is due to the fact that the time of each long-cadence 
measurement was corrected to the Solar System Barycenter. 

As another caveat, we note that many of the contact bina- 
ries exhibit a pseudo-random walk in eclipse phase as well as 
quasi-periodic behavior with typical amplitudes of ^300 sec 
rms (Tran et al. 2013). In addition, the O-C curves for the 
secondary eclipses in these systems are often anti-correlated 
with the primary O-C curve (Tran et al. 2013). The charac- 
teristic timescales for these cyclic changes in phase can range 
from weeks to many months. Therefore, one should be cog- 
nizant of the possibility that O-C periods of the order of the 
3-year Kepler data interval might simply be the lowest promi- 
nent frequency of a random-walk process - especially for con- 
tact binaries. In this work we remain mindful of this possibil- 
ity. We therefore generally require two full orbital cycles (i.e., 
with P tr i p < 600 days) that are strictly periodic before we are 
reasonably confident that a binary is also a good triple-star 
candidate. However, our collection of 39 triple-star candi- 
dates does contain nine systems with P tr j p > 600 days (six of 
these have Pbj n < 1-1 day; three are classified as 'contact bi- 
naries'). The reader can be the judge of the validity of these 
candidates. 

3.3. Candidate Triples 

After eliminating as many false positives as we were able, 
we were left with a list of 39 candidate triple-star systems 
with convincing eclipse timing variations ("ETVs"). The Ke- 
pler Input Catalog (KIC; Batalha et al. 2010) numbers of our 
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Figure 1. O—C data and model fits for 9 systems with KIC numbers between 3228863 and 5376552. The red curves are the total model O—C values. Dark blue 
is the model fit for the Roemer delay (eq. j6|). Light green curves represent the total physical delay (sum of eqs. (8} and J5}). Note that the vertical scales are 
different on all of the plots; the amplitudes of the O — C curves range from a low of 30 sec to a high of 1000 sec. The linear and quadratic terms in the fit have 
been subtracted before the plot is made. 



39 candidate triple stars are summarized in Table 1, along 
with other properties of the targets that are provided in the 
KIC. Among other parameters, we list the orbital period of 
the binary, the Kepler magnitude {K p ) and T e ff of the inte- 
grated light from the system, the depths of the primary and 
secondary eclipses, the mass ratio and "third light" parameter 
(as fou nd w ith the Phoebe binary light curve emulator; see 
section |6.6 1, and an approximate binary orbital eccentricity 



(taken from the Slawson et al. 201 1 catalog). 

The O — C curves for all 39 of the candidate triple systems 
are shown in Figs. [1] [2j [3] HJ and [5 As the reader will see, 
there is a great variety or shapes, of amplitudes, and of sta- 
tistical quality. These, and formal model fits to them, are dis- 
cussed in detail in the following sections. In general, the rms 



deviations from the best fitting curves are in the range of 20 
to 100 sec. The amplitudes of the O — C curves range from a 
minimum of 30 sec to a maximum of nearly 6000 sec. The 
inferred orbital periods of the triples range from 48 days to 
959 days. 

4. SOURCES OF ETV DUE TO THIRD STARS 

4.1. General Expressions 

An eclipsing binary can be thought of as a clock, where the 
clock "ticks" are the binary eclipses. If the binary is circu- 
lar and isolated in space, then the arrival times of the eclipse 
events at the solar system barycenter occur at a constant rate - 
assuming that the binary orbit is neither decaying nor expand- 
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Figure 2. O—C data and model fits for 9 systems with KIC numbers between 5384802 and 7690843. The red curves are the total model O—C values. Dark blue 
is the model fit for the Roemer delay (eq. J6|). Light green curves represent the total physical delay (sum of eqs. (8j and J5}). Note that the vertical scales are 
different on all of the plots; the amplitudes of the O — C curves range from a low of 50 sec to a high of 5000 sec. The linear and quadratic terms in the fit have 
been subtracted before the plot is made. 



ing. When the binary is part of a hierarchical triple system, 
where both the binary and the third star orbit their common 
center of mass, the clock "ticks" are no longer regular. There 
are two basic effects that cause these eclipse arrival times to 
deviate from the pattern of a regular clock, on the timescale 
of the orbital period of the triple. 



4.1.1. Roemer delay 

The first important effect is the classic Roemer delay (or 
light travel time delay) that results from the changing pro- 
jected distance along the line of sight of the center of mass of 
the binary from the center of mass of the triple system. The 
expression for the contribution to the O — C curve from the 



Roemer delay, TZ(t), is 
gg) 

^Roem 



(l- 



? I/ 2 

e ) sinMCOsw + (cosM-e)sincj 



(6) 



where u(t) is the eccentric anomaly, u> the longitude of perias- 
tron, and e the eccentricity, all describing the orbit of the CM 
of the binary about the CM of the triple system. The ampli- 
tude of the Roemer delay is: 



l Roem — 



G l/3 
c(27r) 2 /3 P "'ip 



2/3 



Af 3 sin /trip 
(M 3 +M bin ) 2 /3 



(7) 



where M^ n is the mass of the binary, / tr i p the inclination of the 
orbital plane of the triple system with respect to the plane of 
the sky, and P tr i p is the orbital period of the triple. 
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Figure 3. O—C data and model fits for 9 systems with KIC numbers between 7873703 and 8904448. The red curves are the total model O—C values. Dark blue 
is the model fit for the Roemer delay (eq. j6|). Light green curves represent the total physical delay (sum of eqs. (8l and J5}). Note that the vertical scales are 
different on all of the plots; the amplitudes of the O-C curves range from a low of ~ 150 sec to a high of 6000 sec. The linear and quadratic terms in the fit have 
been subtracted before the plot is made. 



A diagram showing the triple system geometry is given in 
Fig.[6](where some of the quantities labeled appear only in the 
physical delay function - see below for definitions). 

4.1.2. Physical delay 

The second major effect that results in the eclipse timing 
variations is the so-called "physical delay". This results from 
physical changes to the clock, i.e., actual variations in the bi- 
nary period, caused by the third body. Qualitatively, the pres- 
ence of the third body causes the orbital period of the binary to 
be longer than it would be in isolation. The perturbed binary 
period depends on the instantaneous distance from the center 
of mass of the binary to the third star, r tl -i p , and is longest when 
r tr i p is smallest. If the third star is in a circular coplanar orbit, 
the instantaneous distance r tl ; p is a constant, and there are no 



first order effects to be observed in the eclipse times since the 
lengthened binary period is then a constant as well (here we 
are still assuming a circular inner binary orbit). However, if 
the orbit of the third star is either eccentric or inclined with 
respect to the orbital plane of the binary, then the distance 
between it and the binary CM and/or the tidal interaction is 
constantly changing, and so is the binary orbital period. This 
leads to a very distinctive O—C curve. 

A number of approximate analytic expressions have been 
developed for the case of a third body perturbing the orbit 
of a circular binary (see, e.g., Brown 1936; Harrington 1968; 
1969; Soderhjelm 1975, 1982, 1984; Borkovits et al. 2003; 
Agol 2005; Borkovits et al. 201 1) on the timescale of the or- 
bital period of the triple. The perturbative calculation takes 
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Figure 4. O — C data and model fits for 9 systems with KIC numbers between 8938628 and 10613718. The red curves are the total model O — C values. Dark 
blue is the model fit for the Roemer delay (eq. |6j). Light green curves represent the total physical delay (sum of eqs. J8) and J9}). Note that the vertical scales 
are different on all of the plots; the amplitudes of the O — C curves range from a low of 100 sec to a high of 2000 sec. The linear and quadratic terms in the fit 
have been subtracted before the plot is made. 



advantage of the hierarchical nature of the system and ex- 
pands the equations of motion in terms of the small parameter 
£ = TbinAtrip! where rbi n is the instantaneous binary separation 
and r tr i p is the instantaneous distance from the tertiary star to 
the CM of the binary, as defined above. The short period per- 
turbations (those on the timescale of the binary period) are of 
small amplitude (higher order in £) and less interesting ob- 
servationally; averaging over the binary period results in an 
expression for the slower (but higher amplitude) variations in 
the perturbed period of the binary on the timescale of P tr i p . 

The most comprehensive of the expressions for the physical 
delay in the case of circular binaries^ is given in Borkovits et 

9 In this work we utilize two pieces of information to constrain the orbital 
eccentricity of the binaries within our candidate triple stars: (i) analysis of 



al. (2003; but see also Borkovits et al. 201 1 for a more expan- 
sive treatment of perturbations to eccentric binaries). The ex- 
pression there encompasses the perturbations to the period of 
the binary occurring on a timescale equal to P tr ; p , and consists 
of three terms, of which we use two. The two terms appearing 
in the O—C formula which we use are: 



fig) 

^phys 



21- - [<Kt) + e sin <j)(t)- 6(t)] 



(8) 



the epoch-folded light curves (see Table 1 and 8 )6.6) ; (ii) the similarity of the 
O-C curves for the primary and secondary eclipses for the vast majority of 
the systems (especially those with Pb m < 2 days) provides additional evidence 
for the approximate circularity of the binary orbits (see Table 1). 



Rappaport et al. 



XJ 
c 
O 
<j 
m 

o 



O 

I 

o 



— i — i — i — I — i — i — i — i — I — i — i — r 



KIC01 1042923 




0.4 



0.7 



-0.2 



-0.4 




1000 



500 



1000 



Time (days) 



Figure 5. 0-C data and model fits for 3 systems with KIC numbers between 10991989 and 11968490. The red curves are the total model O — C values. Dark 
blue is the model fit for the Roemer delay (eq. J6j). Light green curves represent the total physical delay (sum of eqs. (8) and (9l). Note that the vertical scales 
are different on all of the plots; the amplitudes ofthe O-C curves range from a low of 200 sec to a high of 300 sec. Theunear ana quadratic terms in the fit have 
been subtracted before the plot is made. 

star is too close. In particular, a very close passage of the third 
star could induce a substantial eccentricity in the binary orbit. 
The formulae above, derived assuming a circular binary orbit, 
would then not apply. We find that, for coplanar orbits, the 
formulae agree well with numerical experiments as long as: 



Pi(t) 

^phys 



= (l-Z){sin[20(f)-2v ffl ] 



+ e sin [0(f ) - 2v m ] + - sin [3 0(f) - 2 v m ] } (9) 



where 



iphys ' 



3 M 3 Pl n , ,,-3/2 



3tt M trip P trip 



(10) 



with the following definitions: (f> and 6 are the true and mean 
anomalies of the orbit of the triple system, I is cos 2 i m with 
i m the mutual inclination of the binary orbital plane with re- 
spect to the orbital plane of the triple, and v,„ describes the 
orientation of the periapse of the triple system with respect to 
the binary plane. (See Fig.[6]for definitions of the parameters 
describing the system geometry.) 

The third term in this sequence (not given here), Pj(t), is 
proportional to cot/bi n sin/,„, where ;'bin is the inclination to 
the plane of the sky of the binary orbit. Given that the binaries 
we are studying exhibit eclipses, cot/bin is likely to be small. 
If, in addition, the mutual inclination angle of the two orbital 
planes is small, then the product of cot /bin sin /,„ is likely to 
be negligible for our purposes. Thus, in the present work, we 
exclude this third term. 

As an illustration of how the Roemer and physical delays 
compare, we show in Fig. [7] a plot of the amplitudes of the 
Roemer and physical delays as a function of P tl -i p for six dif- 
ferent assumed periods of the binary. We adopted illustrative 
values of e = 0.3, / tr i p = 60°, and all masses equal to 1 M Q . As 
could be inferred from the analytic expressions, the Roemer 
delay dominates for longer periods of the triple system and 
shorter binary periods, and vice versa for the physical delay. 
The two effects are roughly comparable for a 1-year triple- 
system period and a binary with a 1-2 day period. 

Finally, we note that the accuracy of these analytic ex- 
pressions (eqns. [8] and [9]l has been checked in the original 
Borkovits et al. papers (2003, 2007, 2011) via direct 3-body 
numerical integration. However, one might expect that these 
formulae, derived assuming the parameter £ = flbin Atrip is 
small, must break down if the pericenter passage of the third 
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(1-e) > 5fl b in 



(ID 



Here a tl -i p and e are the full semimajor axis of the triple and 
its corresponding orbital eccentricity, and a\, m is the orbital 
separation of the two stars in the binary. In terms of the orbital 
periods, this corresponds to 



P^-ef' 2 >UP x 



bin 



(12) 



for an assumed set of three equal mass stars. 

An exception to this agreement between the analytic ex- 
pression and the numerical results can occ ur when longer- 
term perturbations (discussed below in §4.2[ ) set in. Since the 
timescales for these longer-term perturbations are typically a 
decade to centuries (see Table 3), they can be fitted (or, effec- 
tively removed) by simply adding linear and quadratic terms 
to the fitting parameters (see 5j5}. 

4.2. Longer-Term Perturbations 

In addition to the perturbations to the orbital period of the 
binary that are discussed above and have a complete cycle 
time equal to the period of the triple system, there are other 
perturbations that occur on typically much longer timescales. 
These include precession of the orbital plane of the binary 
and possible precession of the longitude of periastron of the 
binary, if the binary is eccentric. The approximate timescale 
for these longer-term perturbations is 



Tlongterm C*- 



Pun M 3 



(13) 



(Harrington 1968; 1969; Mazeh & Shaham 1979; Ford, 
Kozinsky, & Rasio 2000; Borkovits et al. 2003; Borkovits et 
al. 2007). Additionally, if the mutual orbital inclination angle 
satisfies 



sin 2 L 



> 2/5 or 39.2° < i m < 140.8° 



(14) 
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Figure 6. Geometry of the triple system. The observer is viewing along the 
+z axis, and the xy plane coincides with the plane of the sky. For the purpose 
of this diagram, as well as for our analysis, we take the binary orbit to be 
circular and its orbital angular momentum vector to lie approximately in the 
xy plane. Of the four angles used in the analysis, i^p, oj, v„,, and i m , the first 

three are indicated in the diagram, while cos im — /*bin '^trip- (Note, however, 
u> = oJiMrd star +•"■■) In words, ( tl ip is the conventional orbital inclination angle 
of the third-star orbit; the mutual inclination angle, !,„, is the angle between 
the two orbital planes; uj is the angle along the outer orbit of the binary CM 
from the plane of the sky to the periastron point; and v„, is the angle along the 
outer orbit from periastron to the plane of the binary. 
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Figure 7. Comparison of the Roemer amplitude (black curve), given by 
eq. l|7) and the physical amplitude (colored curves) given in eq. jlO) as a 
function of the orbital period of the triple system. The various physical de- 
lay curves are for different assumed binary periods, ranging from 0.5 days 
to 20 days, as labeled. See text for a list of the nom inal values that were 
assumed for the other parameters in eqs. l|7) and j I0[ . Dynamically stable 
systems wo uld be expected to lie below ana to the ngh 
(see eq.[l6l. 
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Kozai cycles (Kozai 1962) may set in. In this effect there is 
a cyclic tradeoff between the growth of orbital eccentricity 
of the binary (including when it initially has eu D = 0) and a 
corresponding decrease in i m . If the timescale for this cycle, 
which is the same as Ti ongterm in eq. (13i, is longer than the 
timescales that characterize other perturbations that drive pre- 
cession of the longitude of periastron in the binary, the Kozai 
cycle will not operate (Eggleton & Kiseleva-Eggleton 2001; 
Fabrycky & Tremaine 2007). Moreover, effective damping 
from the two stars in the binary can terminate the Kozai cy- 
cles completely - preferentially leaving i m in the range of 35° 
to 50° (Fabrycky & Tremaine 2007). 

The values of Ti ongterm for all of our triple-star candidates are 
listed in Table 3. They range from ~3 years to 5000 years, but 
with only 7 of the systems having Tiongterm < 15 years. There- 
fore, the generally sinusoidal behavior of these long-term per- 
turbations will look approximately linear or quadratic on the 
3-year timescale of the Kepler data set. And, as a rough ap- 
proximation for representing such behavior, we have included 
a quadratic term in our fit (see ^5}. 

5. ANALYSIS CODE 

5.1. Choice of Fitting Parameters 

Given the above expressions for the Roemer and physical 
delays contributing to the O-C curves, there are a total of 
11 free parameters to fit for, under the assumption that the 
binary orbit is circular. These include 8 parameters which 
describe the triple system as an equivalent binary composed 
of the third star and a star with the mass of the close binary 
at the location of the center of mass (CM) of the close binary, 
and 3 other parameters that describe the O-C curve in the 
absence of the Roemer and physical delays, i.e., a reference 
time, slope, and curvature terms: 

e, eccentricity of the outer orbit, i.e., the triple 
ui, longitude of periastron of the binary CM 
r, time of periastron passage 
i m , mutual inclin. of the orbital planes - eqs 
v m , orientation parameter - eqs. (|SJ, |9]); see Fig 
/trip, orbital period of the triple 



8D, d9 



Mj/Mtnp, mass ratio oc A p h ys (see eq. 10 1 
/(M 3 ) I/3 = cube root of mass function cc A^ oem 
to, reference time 

APbin, mean slope of O-C curve x Pbin 
Pbin, quadratic term 

We have chosen to fit for the mass ratio and cube root of the 
mass function since they are the directly measured quantities 
via the physical and Roemer delays, respectively, if we know 
the orbital period of the triple. The orbital period can gener- 
ally be estimated very well before doing the fit by examining 
the periodicity of the O—C term. The to term is essentially a 
measure of the time of the first eclipse in the sequence. APbin, 
the mean slope of the O-C curve, is not generally zero be- 
cause we used the binary period in the Slawson et al. (2011) 
catalog - based on only 120 days of data - to compute the 
initial set of O—C curves. Finally, the quadratic term could be 
used to measure the orbital decay or expansion of the binary; 
however, we do not expect this effect to be detectable over 
the course of only a few years. Rather, we use this quadratic 
term to take into account possible perturbation s tha t occur on 
timescales substantially longer than P tl -; p (see, ^4.2[ i. 

Depending on the Roemer and physical amplitudes, certain 
among the system parameters may be determined much better 
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than others. For example, if the Roemer delay is dominant and 
the physical delay is negligible, the mass function will be well 
determined but the parameters i m , v m , and M^/M^m will not be 
substantially constrained. On the other hand, if the physical 
delay is well measured but the Roemer amplitude is small, 
then the mass ratio, M^/M^ will be more tightly constrained, 
while the mass function and longitude of periastron will be ill 
defined. 

For a number of reasons we decided against using either 
a conventional Levenberg-Marquardt (LM) or Monte Carlo 
Markov Chain (MCMC) fitting procedure. First, we note that 
there are two different functions (i.e., physical and Roemer 
delays) possibly contributing to the structure of the O-C 
curve, and one does not know, a priori, how much each con- 
tributes. Specifically, in most cases, the two functions are not 
typically orthogonal, and therefore they can trade off against 
one another in the fit. As a result, there can be very large 
regions in parameter space that yield comparably good fits. 
Second, given the large number of systems to deal with, we 
want to search all of parameter space and estimate the uncer- 
tainties at the same time. The LM method is not particularly 
good for exploring parameter space with highly and nonlin- 
early structured correlation functions among the parameters. 
The MCMC fitting techniques is not ideal for exploring wide 
ranges of parameter space, especially when trying to fit 39 
systems. 

We therefore constructed a simpler, though less formal, 
Monte Carlo fitting code that is better suited to the task of 
fitting 39 systems in an automated, hands-off fashion. In this 
approach we choose a random value for each of the follow- 
ing 7 parameters: e, u, r, v m , Pmp, M^/M tT i v , and /(M3) 1 / 3 . 
The parameters are chosen with a uniform distribution over 
their entire plausible ranges. The remaining 4 parameters: i m , 
to, APbin, and Pbin can then be determined via a simple ma- 
trix inversion since they appear linearly in the fitting function. 
(Actually, in the case of i m , it is cos 2 ; m that appears linearly 
in the equations.) 

The uncertainty on the individual data points is determined 
empirically as follows. All data points for a given system are 
assumed to be equally weighted. We then make a first-pass 
run with our simple MC fitting code to find a good set of sys- 
tem parameters. Using that fit, we scale the size of the error 
bars so that the normalized value of the chi-squared statistic, 
x\, is equal to 1. From then on, each time the code is run, we 
use that same value for the error bars on the individual points 
(unless subsequent runs find a substantially improved fit). 

In all subsequent runs, the code operates as follows. If 
the value of xl resulting from a particular selection of pa- 
rameters is xi > 1-3 then we add the ratio of likelihoods, 
ex p[~(X 2 ~~ Xo)/2] (where \o is the value for the best fit), to 
the various probability histograms that are being accumulated 
for each parameter. The code then chooses another random 
set of possible system parameters. If, on the other hand, the 
value of xt resulting from a particular selection of parameters 
is xt < 1-3, then the code does an additional 1000 draws for 
a more restricted range of the parameters surrounding the par- 
ticular choice of parameters that yields the "good" \ 2 value. 
When the 1000 additional draws have been completed, and 
the ratio of likelihoods has been recorded for each draw, the 
broad grid search resumes until another combination of pa- 
rameters is found that yields a value of xt < 1-3- At that 
point, another 1000 localized draws are made, and so forth. 



With this prescription, on average, about half the draws cover 
the broad search while the other half covers a more restricted 
range of parameters. 

This analysis scheme seems reasonably optimum in terms 
of covering all of parameter space while exploring in greater 
detail the regions which yield the best fits. Without full or 
rigorous justification, we also expect it to give approximately 
correct estimates of the parameter uncertainties. 

5.2. The Fitting Runs 

The number of eclipse times, over 13 Kepler quarters, to 
be analyzed in any given binary ranges from only ^40 to as 
many as 2400, depending on the orbital period (except for the 
special case of KIC 10319590 where there are only 19 pri- 
mary eclipses; see Fig. 10 1. The analysis time is essentially 
linearly proportional to the number of eclipses. We chose to 
have the code spend roughly the same amount of time ana- 
lyzing each source rather than drawing the same number of 
random sets of parameters to test. The reason is that for the 
shorter binary periods, the O-C curves become dominated by 
the Roemer delay (since A p h ys oc P^ m whereas AR 0em is inde- 
pendent of fbin)- Since the Roemer delay has one fewer free 
parameter, and is generally simpler in shape than the physical 
delays, such O-C curves can be fit more quickly. 

With this in mind, we typically draw 10 random sets of pa- 
rameters for a fiducial 5-day binary, and this number is scaled 
proportionally to Pbin from that value. The analysis then takes 
a day and a half on a MacBook Air computer, and is adequate 
to yield good fits and system parameters with their uncertain- 
ties. The same analysis was done using 10 6 , 10 7 , and 10 8 
draws (scaled to Pbin/5 days). We found that the 10 7 and 10 8 
draw runs resulted in the substantially the same best fit param- 
eter estimates and any deviations were almost always within 
the 10%-90% uncertainty interval. 

5.3. Test of the Code 

In order to check the basics of the code we simulated eclipse 
timing data for a number of different triple-star systems us- 
ing a 3-body numerical integrator. These include cases where 
the Roemer delay dominated, where the physical delay dom- 
inated, and where the two effects were comparable. White 
noise of rms amplitude equal to 60 sec was added to the sim- 
ulated eclipse arrival times. The artificial data were then an- 
alyzed in exactly the same way as the actual O-C data. The 
results were that the fitting code recovered the correct input 
parameters from the simulation, to within the 10% - 90% er- 
ror constraints (the same as we list in Tables 2 and 3). 

6. RESULTS 
6.1. Overview 

The results of the automated fits to the 39 triple-star candi- 
dates are shown in five multi-panel figures (Figs. [TIB}. They 
are arranged simply in order of their KIC number. In each 
panel, the red curve is the overall fit to the O-C curve, and is 
the sum of the Roemer and physical delays, which are shown 
separately as the blue and green curves, respectively. 

The fitted parameters and their uncertainties are listed in Ta- 
bles 2 and 3 along with the 10% and 90% (upper and lower) 
confidence limits. Table 2 gives, in addition to the binary pe- 
riod, four quantities related to the masses which are derived 
entirely from fitting the O-C curves. These are the mass ra- 
tio, M 3 /M tr i p , the mass function, M\ sin 3 //(M 3 +M b i n ) 2 , and 
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the quantities M3 sin i and Mb; n sin i, derived from the mass 
ratio and mass function - to the extent allowed by the uncer- 
tainties. We also list the amplitudes of the Roemer and phys- 
ical delays (the 10% and 90% probability limits are given in 
curly brackets). 

In Table 3 the remainder of the fitted parameters, eccen- 
tricity, e, and time of periastron passage, r (relevant to both 
Roemer and physical delays), the longitude of periastron, oj 
(appearing in the Roemer delay only), and the mutual orbital 
inclination angle, i m , and orientation angle, v m (both related to 
the description of the physical delay), are given. Table 3 also 
lists the rms of the residuals with respect to the best fitting 
O-C curve, as well as the calculated timescale for longer- 
term perturbations (see eq.fT3i. 

EE 



as well as Table 2 shows that 19 of 



A perusal of Figs, 
the O-C curves are dominated by the Roemer delay, 1 1 are 
dominated by the physical delay, while the remaining 9 ob- 
jects have more competitive Roemer and physical amplitudes 
(here "dominant" is defined as a > 3 : 1 ratio). If "dominant" 
is defined by a ratio of > 5 : 1, then the corresponding num- 
bers are 18 Roemer, 8 physical, and 13 comparable. The Roe- 
mer delay dominated systems all have binary periods of < 2 
days, consistent with the diagram in Fig. [7] Conversely, all 
the systems with the longer orbital periods (e.g., > 5 days) 
are dominated by physical delays. 

6.2. System Parameter Constraints 

A review of Table 2 will show that for systems that are 
dominated by the Roemer delay, the cube root of the mass 
function is indeed determined with greater fractional accu- 
racy (^10%) than is the mass ratio (typically >40%). This 
follows from the fact that the Roemer amplitude is directly 
proportional to the cube root of the mass function. Addition- 
ally, in this circumstance, the parameters uj, r, and e are all 
relatively well determined, but the parameters strictly associ- 
ated with the physical delay, v m and i m , are generally poorly 
constrained. Conversely, for the systems where the physical 
delay dominates, the mass ratio, Mj/M^p, is determined to 
a substantially better fractional accuracy (~30%) than is the 
cube root of the mass function (typically >50%). Again, this 
is due to the fact that the physical amplitude is directly pro- 
portional to the mass ratio. As well, the parameters i m , r, and 
e, are better determined than oj which is only relevant to the 
Roemer delay. The parameter v m , generally seems not well 
constrained, except in six systems - all ones with dominating 
physical delays. 

One might guess that for those 9-13 systems where the Roe- 
mer and physical delays are more comparable (smaller than 
3:1 or 5:1 ratios, respectively) both the mass ratio and mass 
function could be well determined. This does not appear to 
be the case in practice. The reason is due to the fact that the 
two sets of functions representing these delays are not sub- 
stantially orthogonal, and therefore the two functions can add 
in different ways, consistent with the constraints of the param- 
eters r, u, and v m to produce the total observed amplitude. It 
turns out that the Roemer and physical delays, when compa- 
rable, can vary together in amplitude over a fairly wide range 
while the longitude of periastron, ui, in turn, changes their rel- 
ative phase in such a way that the sum of the two functions 
adds to be roughly a co nstan t (and thereby matches the ob- 
served O-C curve; see §6.3| for details). Thus, in no specific 
system, do we obtain very tight constraints on both M3 sin 3 i 
and Mbin sin 3 i (i.e., with both being determined to better than, 
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Figure 8. Example of the correlation between the eccentricity of the outer 
orbit, i.e., the triple, and the mass ratio, M^/M u \ v for a system in which the 
physical delay dominates: KIC 9714358. The colors are scaled according to 
the relative probability with white and red the highest, blue and purple the 
lowest. 
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Figure 9. Example of the correlation between the cube root of the mass 
function, /(M3) 1 / 3 , and the longitude of periastron, ui, of the outer orbit, 
i.e., the triple, for a system in which the physical and Roemer delays are 
comparable: KIC 9451096. The colors are scaled according to the relative 
probability with white and red the highest, blue and purple the lowest. 



e.g., 20%). 

This kind of correlated behavior between the physical and 
Roemer contributions is no t nearly so strong when only one of 
them dominates (see §6.3| l. The reason is that for given binary 
and triple star periods, as well as eccentricity, the physical 
delay has a tight upper limit that is proportional to M^/M^ 
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which, by definition, can never exceed unity. Since the phys- 
ical delay amplitude is proportional to P£- , while the Roe- 
mer amplitude is independent of Pbin, for short period binaries 
(i.e., < 0.7 days) it becomes difficult for the physical ampli- 
tude to contribute much to the O-C curves, notwithstanding 
any issues of orthogonality. Conversely, the Roemer delay is 
proportional to the cube root of the mass function which, has 
an upper limit of M3. While in principle, it is possible for the 
mass of the third star to take on any value, unless it is a fairly 
evolved giant, it is unlikely to have a mass greater than a few 
Mq since both stellar radius and T e s were constrained by the 
nature of the stars selected for inclusion in the KIC (Batalha 
et al. 2010). Therefore, for systems with long binary orbital 
periods it is difficult for the amplitude of the Roemer delay 
to contribute much regardless of how well the shape might 
match the O-C curve. 



6.3. Correlations Among the Parameters 

We have tried to select a convenient, consistent set of sys- 
tem parameters to fit for all of our candidate triple stars, re- 
gardless of whether they are dominated by the Roemer or 
physical delays. It is somewhat inevitable that some of the 
paramet ers c an become substantially correlated (see discus- 
sion in 5 6.2 1 when the physical delay dominates, vice versa, 
or even when the two effects are comparable. Here we show 
two examples of this type of correlation taken from our Monte 
Carlo fitting code. In Fig.lSlwe show the correlation between 
the eccentricity of the triple (i.e., the outer orbit) and the mass 
ratio, M 3 /M trip , for the example of KIC 9714358 which is 
dominated by the physical delay. In the case of physical delay 
only, the amplitude is roughly proportional to the product of 
these two quantities, and we then expect just such a correla- 
tion as is seen in Fig.jS] This can be shown analytically for the 
case of coplanar orbits from eq. ([8]) where the term in square 
brackets on the right hand side, [0(f) + e sin 4>(t)-9(t)], can 
be expanded in a series for small eccentricities as ^3e sin <fi(t) 
(Murray & Dermott 2000), while the M 3 /M ai? part of the pro- 
portionality is found in eq. ( ff0| . For non-coplanar orbits, one 
of the terms in eq. Q is not proportional to e while the other 
two terms are; therefore, the correlation becomes less pro- 
nounced as the mutual orbital inclination increases. 

We now consider the key correlation for the case where 
the physical and Roemer delays are more comparable. In 
Fig. [9] we show the correlation between the cube root of 
the mass function, /(M3) 1 / 3 and the longitude of periastron, 
lj, for the case of KIC 9451096. The correlation seen in 
Fig. [9] is quite strong and symmetric around 180°. The zero 
delay point of the physical delay typically occurs near the 
time of periastron passage, r (especially as i m — > 0), while 
the Roemer delay is zero at ~ t - loP^/Iit. Therefore, if 



1.02 



/i Roem ~ /iphy S the two functions will have a combined ampli- 
tude A -c — 2A Roem | cos(w/2)|. It then follows that A Roem ~ 
Aphys — jAo-c/ 1 cos(oj/2) I and these two parameters (AR 0em 
and id) are thus highly correlated, as seen in Fig. [9] 

6.4. Dynamical Stability of Orbits 

We mention in passing that, as a sanity check on the or- 
bital solutions we have found, the mutual orbits of the three 
stars would be expected to have long-term dynamical stabil- 
ity. The stability criteria for triple systems have been studied 
for decades, and are conveniently summarized by Mikkola 
(2008). In particular, we cite here the expression due to 
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Figure 10. Example of a system (KIC 10319590) where the eclipse depths 
exhibit strong variations with time. In this extreme case, the eclipses com- 
pletely disappear after ~400 days, presumably due to the precession of the 
binary orbital plane caused by the presence of the inferred third body. 



Mardling & Aarseth (2001): 
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where, again, e is the eccentricity of the triple system. Ex- 
pressed in terms of the orbital periods, this stability criterion 
comes to: 

1/10 (l +e ) 3 /5 



Ari P >4.7 



M, 



trip 

M bin 



(l-e) 9/5 



(16) 



Note that, while we do not know the masses of the binary and 
triple very accurately, the dependence on masses in eq. ( [T6| 
is extremely weak. Moreover, in most cases we have a good 
handle on e, and an excellent measurement of both Pbin and 
P tl i p , so we can then directly compute that all of our triple- 
star candidates are nominally stable. This is another sanity 
check that suggests that these are true triple stars and not false 
positives, since false positives should not be biased towards 
satisfying stability requirements. 

6.5. Supplemental Information Required 

The triple-star candidates identified in this work from the 
Kepler data base will require supplemental information in or- 
der to reasonably infer their full sets of system parameters 
with astrophysically useful accuracy. For some of the sys- 
tems there can be up to three pieces of supplemental informa- 
tion from the Kepler light curves themselves. It is beyond the 
scope of this paper to try to utilize this information, but we list 
them here for the interested reader. Seven of the systems ex- 
hibit secularly varying eclipse depths (see Table 1). The most 
extreme case of secularly varying eclipse depths is the case of 
KIC 10319590 whose flux vs. time is shown in Fig. 10 where 
the eclipses disappear after ^400 days. Two of the systems 
show eclipses of, and/or by, the third body (Carter et al. 2013). 
Finally, four of the systems have O- C curves for the primary 
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Figure 11. Example of a system (KIC 7955301) where the O — C curves for 
the primary and secondary eclipses lie on a "divergent" path - at least for the 
3-year duration of the Q1-Q13 data. As well, the two O—C curves even have 
somewhat different profiles. 

and secondary eclipses that are different in shape and/or sys- 
tematically shift in phase with respect to one another. A good 
example of the later effect is exhibited in Fig.fTTlfor the case 
of KIC 7955301 where the O — C curves for both the primary 
and secondary eclipses are shown. In total, however, there are 
only seven systems, of the 39, which exhibit one or more of 
these three different features. (See Table 1 for a summary.) 

In these seven cases, the supplementary information from 
the Kepler photometry can be modeled with a 3-body code 
to gain a much more complete understanding of the system 
parameters (see, e.g., Carter et al. 201 1 and 2013). 

For these seven systems, as well as the remainder of the 39 
triple-star candidates, it will be important to obtain Doppler 
spectroscopy. Even a high-quality, single-epoch spectrum, 
could provide significant insight into the nature of the three 
constituent stars. Measuring the Doppler motion of the bi- 
nary, and, even better, all three stars, would lock in most of the 
physically important system parameters that are only loosely 
constrained through the analysis of the O—C curves alone. 

In general, the binary orbital periods are quite short (only 
seven have P D in ^ one week), so it will not take much ob- 
serving time to unravel the properties of the binary (e.g., its 
masses and luminosity contribution to the triple). The periods 
of most of the triples range from 48 days to 1 year. The me- 
dian period is ^330 days. Therefore, Doppler spectroscopy 
aimed at determining the properties of the triple orbit would 
have to span a good portion of the observing season for the 
Kepler field. 

6.6. Binary System Light Curves 

To gain some further insight into the constituent stars in the 
39 systems we have identified, we have constructed folded 
light curves for each of the binary stars in these systems. We 
then used the Phoebe binary light curve modeling code (Prsa 
& Zwitter 2005) to fit the binary system parameters, allowing 
for the "third light" parameter (presumably largely due to the 
light contribution of the third star) to be a variable. The results 



for both the contribution of the "third light" and the mass ratio 
of the two stars in the binary, qu n , are listed in Table 1. In 
principle, this information can be used in conjunction with the 
constraints on M3 and M avp found from the analysis of the O— 
C curves (see Table 2) to infer the three masses individually, 
albeit with wide uncertainties. 

We were also able to use the Phoebe fits to check the orbital 
eccentricities of the binary systems as reported by Slawson et 
al. (2012), and we find reasonable agreement, though with the 
Phoebe values of e D in tending to be a bit lower. The value 
of e D i n is important for the expected form of the physical de- 
lay curve; the O — C curves can be noticeably affected when 
e D in >0.05 or so. Table 1 lists the binary eccentricities com- 
puted from values given in the Slawson et al. (2012) catalog, 
but replaced in four cases with the Phoebe result (where the 
former value of e D in was more than 3 times higher than the 
Phoebe value). In all, six of the systems have e D m It 0.075, 
and we note that the fitted triple-star parameter values for 
these could be significantly different from the true system pa- 
rameters. 

7. DISCUSSION 

In all, we computed and examined the O- C curves for some 
2000 Kepler binaries. We found that approximately 50% of 
these yielded quite useful portraits of the source eclipse tim- 
ing behavior, with typical rms scatter of better than 100 sec- 
onds. Some 20% were contact (or otherwise short-period) bi- 
naries that tended to exhibit erratic, or random-walk like be- 
havior that made it difficult to search for periodic signatures 
of third bodies. The remaining 30% yielded minimal or no 
useful information on the eclipse timing. In some cases this 
latter category could be attributed to eclipse depths that were 
too small, stellar noise (i.e., starspots, stellar oscillations, etc.) 
that was not sufficiently filtered out, and/or inadequacies in 
our eclipse detection algorithm^] We believe that the 50% of 
binaries for which we were able to obtain good eclipse timing 
information is sufficient so that our findings are not substan- 
tially biased. 

Notwithstanding the above general statements about our 
search, there are quite a few observational selection effects in 
play. These include the construction of the Kepler input cat- 
alog itself which selected for certain spectral types and radii. 
Then, there is the binary detection efficiency for the various 
stars within the KIC. Among other things, this depends on 
stellar pulsations and starspot activity. Within our search for 
triples, the depth of the eclipses, which in part depends on the 
brightness of the third star, affects the timing accuracy. The 
erratic timing behavior of many contact binaries (at the ^300 
sec rms level) makes it harder to detect tertiary companions 
(via eclipse timing variations) in these systems. Finally, if we 
limit ourselves to seeing 1.5-2 orbital cycles of the triple, 
then orbital periods greater than ^900 days are nearly ruled 
out. In fact, in our visual inspection of the set of O — C curves 
we see numerous such potential longer-period triple-star can- 
didates (see also Gies et al. 2012). On the short period end, 
there are many beat periods, between the Kepler cadence and 
the binary period, up to ^20 or 30 days, thus, a number of 
shorter triple-star periods would be hard to have confidence 
in. 

The periods of the triple-star candidates we found are plot- 

10 The fraction of systems (~30%) that yielded no useful O-C curves 
did not improve with the use of our newly developed, more formal cross 
correlation analysis (mentioned earlier in the text.) 



14 



Rappaport et al. 



1000 



3h 



100 



T3 

o 




Binary Period (days) 

Figure 12. Plot of the periods of the candidate triple systems vs. the period 
of the binary system they contain. The blue line indicates the locus of points 
where Ptrip/fljin = 10, as a representative stability criteri on. Most systems 
should lie to the left of this line which is taken from eq. ijjjl with e = 0.3 



The horizontal green line is a rough lower limit to values of P tl j p that can 
be detected via the Roemer delay with the Kepler Q1-Q13 data set, given a 
sensitivity of ~50 seconds in detectable amplitude (see eq. 17). Finally, the 
red line is a rough upper limit to values of P t rip that can be detected via the 
physical delay with the Kepler Q1-Q13 data set given a sensitivity of ~50 
seconds in detectable amplitude (see eq . 1 1 0) . An assumed value of e = 0.3 
was used to evaluate this latter limit. Systems to the left of the red line are 
typically detected via the Roemer delay. 



ted vs. their binary periods in Fig.[T2] We show a rough dy- 
namical stability bound on the rightlolue curve). This limit 
is derived from eq. ( 16 1 for an assumed typical eccentricity of 
the triple equal to 0.37Most of the triples should lie to the left 
of this curve. If we assume a typical sensitivity in the O-C 
curves of ~50 sec, the corresponding triple period required 
to produce a detectable signal purely via the Roemer delay is 
about 45 days (see green curve in Fig. [12] above which we 
should be able to detect the light-traveRime effects). Here 
we have assumed all 1 M constituent stars, and an orbital 
inclination of the triple equal to 60° (see eq.[7]i. The limiting 
triple-star periods for the physical delay are indicated crudely 
by the red curve in Fig.[T2]>. This is based on eq. ( [T0| with 
e = 0.3 and all equal constituent masses. Systems detected via 
the physical delay should lie to the right of this curve for an 
amplitude sensitivity in the O-C curves of ^50 sec; systems 
to left of this line are detected via the Roemer delay. Finally, 
it is difficult at best to confirm any triples with P tl ; p > 1000 
days (see also Gies et al. 2012). 

Thus, Fig.[T2]indicates that most of the 39 triple-star candi- 
dates are reasonably well dispersed (in log space) around the 
zone of detectability and stability. 

Because of the various observational and analysis selec- 
tion effects alluded to above, it is difficult for us to draw 
far-reaching conclusions about the fraction of binary systems 
with relatively close tertiary companions. However, there are 
some things we can say in this regard. There were approxi- 
mately 1000 of the Kepler binaries which yielded useful con- 
straints on the eclipse timing via our particular approach to 
the analysis. There were some 39 triple star candidates found 
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Figure 13. Distribution of the mutual inclination angles of the 39 triple-star 
candidates. Note that, in general, the uncertainties in i m are larger than the 5° 
bin size used for the histogram. 



among these with 48 < Ptrip 900 days, spread roughly uni- 
formly in logP tr ip. Without trying to be too precise, we can 
say that we see evidence for roughly a comparable number 
of potential candidate triple systems with P tl i p in the range of 
^1000-2500 days, where only at most one to a fraction of 
an orbital cycle is revealed. This would suggest that perhaps 
^8% of close binaries have tertiary companions that have or- 
bital periods of less than ~7 years. Again, the O-C sensitivity 
limit here is ^50 sec (rms scatter) with which we are able to 
time the eclipses. 

Finally, in terms of the completeness of our initial survey 
for triple systems, we note that some of the companions to 
binaries with P\, m < 1 day and P tr j p < 30 days can produce de- 
lays that are too small (i.e., less than a few tens of seconds) to 
be detectable with the current approach. In particular, note the 
unpopulated region in the bottom lower left corner of Fig. 12 

Among the most popular formation theories for very close 
binaries (e.g., with Pbin ;$ 3 days) are those which invoke a 
third star, even if quite distant (with P tl i p up to 10 yr), to 
effect the closeness of short-period binaries. These scenar- 
ios typically involve so called "KCTF" (Kozai cycles with 
tidal friction; Eggleton & Kiseleva-Eggleton 2001; Fabry cky 
& Tremaine 2007; but it is also possible that magnetic brak- 
ing plays a role, e.g., Verbunt & Zwaan 1981; Matt & Pudritz 
2005). Fig. [13] shows the distribution of mutual orbital in- 
clination angles in our sample of triple-star candidates. This 
distribution was produced without regard for the large uncer- 
tainties in the measurements of i m which typically exceed the 
bin width of 5° used here. Nonetheless, there is something of 
a very suggestive peak in the mutual orbital inclination range 
of 35° -45° predicted by Fabrycky & Tremaine (2007) for the 
KCTF scenario. Within our parameter uncertainties, it is quite 
possible that the Kozai cycle is no longer operative in any of 
these systems. 

The present study of tertiary stars around short period bi- 
naries is quite complementary to those of others (see e.g., 
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Mazeh 1990; Tokovinin et al. 2006; Pribulla & Rucinski 
2006; D'Angelo, van Kerkwijk, & Rucinski 2006; Rucin- 
ski, Pribulla, & van Kerkwijk 2007; Raghavan et al. 2010). 
In particular, the region of orbital period space covered by 
Tokovinin et al. (2006; 20 < P trip < 10 5 yr and P hia < 25 days) 
is almost exactly complementary to ours which extends up to 
P U i p < 3 yr and covers the same range of binary periods (see 
Fig. 12 in Tokovinin et al. 2006). If we somewhat arbitrar- 
ily adopt a distribution of triple periods that is constant per 
logarithmic interval, then our detection of ^4% triples over 
a factor of 20 in P tr i p (1.3 dex) is consistent with a signifi- 
cant fraction of all close binaries having tertiary companions 
(Tokovinin et al. 2006; Pribulla & Rucinski 2006; Raghavan 
et al. 2010). If we assume that possible triple-star periods 
cover ~20 days - 10 5 yr (6.3 dex), then we have examined 
~l/5 of this range. Therefore, we might speculatively ex- 
trapolate our results to suggest that >20% of close binaries 
have tertiary companions. Tokovinin et al. (2006) find a much 
higher fraction for binaries with Pbin ^ 3 days, and a more 
comparable one to our value for fbin ^12 days. Thus, given 
all the uncertainties, our results may not be dissimilar. How- 
ever, we do not have the statistics to comment on the tertiary 
fraction separately for binary periods above and below this 
transition period of ^10 days (see, in particular, Fig. 14 of 
Tokovinin et al. 2006). 

8. SUMMARY AND CONCLUSIONS 

We have analyzed the Kepler binary data set for eclipse tim- 
ing variations, with the intention of identifying the presence of 
third bodies. We found some 39 plausible candidates for triple 
star systems, eight of which had been previously found by the 
members of the Kepler team, but only a few of these had been 
studied in any detail. Some were found via tertiary eclipses, 
while others were detected from systematic variations in their 
O-C curves (in the latter case using typically only ~ 1/10 of 
the data in the current study). We have subjected all of the 39 
systems in this study to an analysis which includes possible 
Roemer delays as well as physical delays. All the best fits are 
physically sensible, if not necessarily the final solutions. 

We have shown that at least 8% of close binaries have ter- 
tiary companions with P t rip ^ 7 years. This is in agreement 
with other surveys covering tertiaries in much wider orbits 
over a larger dynamic range in periods. 

In order to fully determine the system parameters in 
the triple-system candidates we have found, ground-based 
Doppler spectroscopy will be required. This is already be- 
ing pursued for a number of the systems (see, e.g., Carter et 
al. 2011; 2013). Moreover, for those systems which exhibit 
other effects of the third body, such as tertiary eclipses, vary- 
ing binary eclipse depths, and/or the effects of binary eccen- 
tricity, there is need for analysis with a 3-body dynamics code. 
We consider our list of triple-star candidates something of a 
starting point for such more extensive studies, both observa- 
tionally and in modeling. 

We were gratified to find that this exercise has proven a very 
good way of finding non-eclipsing triples. 
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Table 1 

Candidate Triple-Star Systems Found in Kepler Data Base 



Source K p l T^f' Prim. Eel. Sec. Eel. et,in 3 <?bin 4 ^/'-trip 4 Vary. Eel. Tertiary Diverg. Prim. 
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3228863 s 


0.730942 


11.82 


6561 


0.440 


0.220 


0.034 


1.20(1) 


_ 


_ 


_ 


no 


4647652 


1.064820 


10.70 


6265 


0.077 


0.021 


0.078 


0.24(1) 


0.224(4) 


_ 


_ 


no 


4909707 


2.302370 


9.20 


NA 


0.043 


0.018 


0.073 


0.075(1) 


0.163(3) 


_ 




_ 


4940201 


8.81659 


13.16 


5284 


0.027 


0.013 


0.083 


0.045(1) 


0.189(1) 


_ 


_ 


no 


5039441 


2.151390 


11.50 


5943 


0.259 


0.019 


0.036 


0.72(1) 


0.018(2) 


_ 


_ 


no 


5128972 


0.505317 


13.23 


5776 


0.094 


0.047 




0.53(2) 


0.207(2) 


_ 


_ 


no 


5264818 


1.905052 


8.70 


9212 


0.013 


0.011 


_ 


1.43(1) 


_ 


_ 
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6520 
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0.109 


_ 
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_ 


_ 


no 


5384802 
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13.70 


6433 


0.020 


0.020 


0.072 


0.42(1) 


0.076(5) 


_ 


_ 


no 


5771589 


10.74007 


11.81 


5927 


0.0011 


0.0007 


0.010 7 


0.03(1) 


0.013(1) 


ves 




ves 


6370665 


0.932316 


14.00 


7386 


0.090 


0.075 




0.52(1) 


0.081(32) 






no 


6525196 


3.42060 


10.15 


5966 


0.162 


0.147 


0.038 


0.71(1) 


0.024(1) 






no 


6531485 


0.676991 


15.55 


5587 


0.021 


0.017 


0.048 


0.032(1) 


0.084(1) 






no 


6545018 


3.99146 


13.75 


5594 


0.291 


0.226 


0.075 


0.77(1) 








no 


7289157 


5.26640 


12.95 


5922 


0.062 


0.006 


0.064 


0.10(1) 


0.299(1) 


yes 


yes 


yes 


7668648 


27.8184 


15.32 


5875 


0.232 


0.094 


0.074 


0.49(1) 


0.014(2) 


yes 


yes 


yes 


7690843 


0.786259 


11.08 


4827 


0.049 


0.020 


0.059 


0.05(1) 


0.303(1) 


_ 




no 


7837302 


23.83530 


13.7 


NA 


0.026 


none 


0.17 


0.010(1) 


- 


- 


- 


NA 


7955301 


15.3266 


10.38 


4821 


0.016 


0.01 


0.20 


0.23(1) 


0.031(2) 


yes 


- 


yes 


8023317 


16.57828 


12.9 


5625 


0.034 


0.002 


0.057 


0.15(1) 


< 0.001 


yes 


- 


- 


8043961 


1.559210 


10.74 


6348 


0.207 


0.170 


0.028 


0.62(1) 


0.140(1) 


- 


- 


no 


8192840 


0.433547 


13.47 


6136 


0.033 


0.028 


- 


0.61(1) 


0.279(3) 


- 


- 


no 


8386865 


1.25800 


12.0 


8510 


0.005 


0.005 


0.59 


0.053(3) 


- 


- 


- 


- 


8394040 


0.392128 


14.46 


5697 


0.042 


0.034 


- 


1.15(2) 


0.53(1) 


- 


- 


no 


8719897 


3.15142 


12.39 


4906 


0.195 


0.176 


0.061 


0.23(1) 


0.015(3) 


- 


- 


no 


8904448 


0.865981 


13.88 


7820 


0.180 


0.049 




0.31(1) 


0.065(6) 








8938628 


6.86219 


13.68 


5602 


0.050 


0.034 


0.062 


1.42(1) 


0.037(1) 


yes 




no 


9451096 


1.25039 


12.64 


NA 


0.233 


0.087 


0.063 


0.46(1) 


0.062(1) 






no 


9714358 


6.47418 


15.00 


4825 


0.185 


0.012 


0.041 7 


0.36(1) 


0.031(1) 








9722737 


0.418528 


14.93 


6517 


0.102 


0.088 




0.50(1) 


0.119(4) 






no 


9912977 


0.943916 


13.7 


NA 


0.292 


< 0.015 


0.01 7 


0.20(1) 










10095512 


6.01720 


13.05 


5795 


0.113 


0.051 


0.082 


0.77(1) 


0.030(1) 






no 


10226388 


0.660658 


9.98 


NA 


0.174 


0.131 




0.18(1) 








no 


10319590 


21.3216 


12.13 


5518 


0.026 


0.008 


0.108 


0.40(1) 


0.079(1) 


yes 






10613718 


1.175880 


10.85 


5080 


0.006 


0.005 


0.099 


0.05(1) 


0.016(1) 






no 


10991989 


0.974475 


10.28 


5021 


0.008 


0.004 


0.05 7 


0.007(1) 


0.167(1) 






no 


1 1042923 


0.390164 


14.32 


6086 


0.210 


0.208 




0.48(1) 


0.153(2) 






no 


1 1968490 


1.078899 


13.70 


NA 


0.033 


0.017 


0.052 


0.043(1) 


0.228(1) 






no 



Note. — (1) The Kepler magnitude and effective temperature are taken from the Kepler input catalog; (2) Depths of the primary and secondary eclipses, based on our 
epoch-folded light curves; (3) Eccentricity of the binary, taken from Slawson et al. (2011) as e b in = [(esinajbin) 2 + (ecosw b in) 2 ] l ' 2 , except where otherwise noted; (4) 
Mass ratio of the binary stars, q bi „, and the fraction of the total Kepler luminosity contributed by the third star, Li/L,n f , as analyzed with the Phoebe binary light curve 
fitting code (the number in parentheses reflects the statistical uncertainty in the last significant digit(s)); (5) See Table 3 for references; (6) This object is the same as the 
eclipsing binary V404 Lyr (see, e.g., Pigulski et al. 2009); (7) Substituted with values from our Phoebe light curve analysis. 
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Table 2 

Candidate Triple-Star Systems Found in Kepler Data Base 



Source 




^trip 


M 3 /M tnp 


f(M 3 ) 2 


M 3 sin 3 imp 


M bin sin 3 / tti p 


ARoem"^ 


Aphys 4 




(days) 


(days) 






(M Q ) 


(Mq) 


(Mq) 


(sec) 


(sec) 


3228863 


0.730942 


668.4 


42 




0.017(0.016,0.017} 


0.10(0.07,0.28} 


0.13 


0.08,0.90} 


189(187,194} 


3.5(2.0,4.0} 


4647652 


1.064820 


753.5 


41 


U.ZO,U.J 3 j 


0.023(0.012,0.039} 


0.13(0.08,0.31} 


0.17 


0.09,0.80} 


228(183,274} 


7.5(4.7,10.4} 


4909707 


2.302370 


505.3 


70 


n n S6 i 

U. JU,U.oO J 


0.510(0.230,1.053} 


1.08(0.47,2.65} 


0.40 


0.11,2.18} 


493(378,627} 


122(81,189} 


4940201 


8.81659 


361.6 


52 


35 77 } 


0.268(0.042,1.266} 


1.08(0.19,3.22} 


0.80 


0.14,3.33} 


318(171,534} 


1209(846,1768} 


5039441 


2.151390 


667.8 


42 


n 76 n ^7 1 

U.ZO,U.J / ) 


0.026(0.011,0.061} 


0.15(0.08,0.36} 


0.17 


0.09,0.81} 


220(163,293} 


39(24,60} 


5128972 


0.505317 


447.8 


U.J J 


n m n 6Q i 


0.094(0.079,0.108} 


0.29(0.20,0.66} 


0.23 


0.09,1.08} 


259(244,271} 


3.9(2.7,4.9} 


5264818 


1.905052 


296.3 


42 


n 76 n 6n\ 


0.037(0.015,0.094} 


0.21(0.09,0.66} 


0.24 


0.09,1.69} 


145(107,196} 


66(42,99} 


5310387 


0.441669 


214.2 


16 


n inn mi 

U. 1 U,U.ZU [ 


< 0.001 


0.03(0.02,0.07} 


0.15 


0.09,0.55} 


3 1 ( 27, 37 } 


2.4(1.5,3.7} 


5376552 


0.503819 


334.5 


32 


U.ZU,U. 3y j 


0.008(0.007,0.009} 


0.08(0.06,0.19} 


0.16 


0.09,0.72} 


94( 91, 98} 


3.3(2.0,4.1} 


5384802 


6.08309 


254.8 


41 


27 71} 


0.075(0.007,0.972} 


0.48(0.07,2.68} 


0.53 


0.11,2.68} 


165( 75,387} 


754(559,1168} 


5771589 


10.74007 


113.2 


n ^5 


U.JZ,U.Jo ) 


0.073(0.009,0.247} 


0.59(0.08,2.08} 


1.10 


0.15,3.96} 


95( 48,142} 


4193(3913,4493} 


6370665 


0.932316 


285.9 


76 


U. 1 / ,U. 3 A \ 


0.004(0.003,0.005} 


0.06(0.04,0.15} 


0.15 


0.08,0.72} 


67{ 61, 74} 


9.0(5.7,10.9} 


6525196 


3.42060 


415.8 


38 


n 07 n i 

U.Z/.U. Jo) 


0.063(0.031,0.201} 


0.59(0.15,1.33} 


0.85 


0.16,3.45} 


215(171,318} 


127(91,189} 


6531485 


0.676991 


48.3 


61 


n 34 n 77 1 

U.Jt,U. / / j 


0.173(0.014,0.613} 


0.32(0.13,2.77} 


0.18 


0.08,3.22} 


72( 31,109} 


83(58,109} 


6545018 


3.99146 


90.6 


29 


u.zu,u.tu j 


0.038(0.005,0.297} 


0.51(0.07,1.87} 


1.21 


0.16,3.78} 


66( 33,131} 


572(439,866} 


7289157 


5.26640 


243.8 


52 


n n 76 \ 

U. JU.U. /o ) 


0.187(0.021,1.065} 


0.74(0.14,2.83} 


0.57 


0.12,2.90} 


218(104,387} 


737(504,1029} 


7668648 


27.8184 


203.7 


10 


U.Uo,U. 1Z J 


0.001 {< 0.001,0.004} 


0.14(0.02,0.41} 


1.31 


0.22,3.65} 


37{ 21,55} 


4759(4097,5401} 


7690843 


0.786259 


74.3 


40 


06 64! 


0.071(0.026,0.147} 


0.41(0.13,1.05} 


0.49 


0.11,2.80} 


71{ 51, 91} 


40(24,61} 


7837302 


23.83530 


959.3 


44 


U.ZO,U. / J ) 


0.177(0.017,1.281} 


1.03(0.15,3.37} 


1.13 


0.16,3.74} 


528(244,999} 


2770(1748,4545} 


7955301 


15.3266 


209.5 


U.JO 


fin film 

U. JZ,U. JV J 


0.094(0.012,0.277} 


0.73(0.10,2.18} 


1.30 


0.18,3.96} 


156( 79,223} 


5788(5464,6131} 


8023317 


16.57828 


613.5 


10 


n Asi n i zli 

U.Uo,U. 1^/ 


0.001 (< 0.001,0.007} 


0.10(0.02,0.42} 


0.85 


0.17,3.33} 


70( 41,131} 


528( 410, 680} 


8043961 


1.559210 


476.7 


41 


U.ZJ.U. JO ) 


0.034(0.028,0.045} 


0.21(0.12,0.49} 


0.29 


0.10,1.42} 


194(179,213} 


24(15,33} 


8192840 


0.433547 


803.9 


38 


n n 47 1 

U.Xj,U.*W ) 


0.015(0.011,0.019} 


0.10(0.07,0.26} 


0.16 


0.09,0.85} 


208(187,223} 


1.9(1.3,3.1} 


8386865 


1.25800 


293.0 


n 55 


n 36 n 67 1 

U.JO.U.O / ) 


0.063(0.047,0.117} 


0.23(0.14,0.62} 


0.18 


0.08,1.08} 


171(156,210} 


37(26,49} 


8394040 


0.392128 


394.8 


U. / 1 


nil fi fiii 


0.353(0.287,0.414} 


0.70(0.50,1.58} 


0.28 


0.10,1.81} 


369(345,391} 


5.4(3.5,7.7} 


8719897 


3.15142 


332.7 


0.52 


0.36,0.70} 


0.158(0.086,0.283} 


0.59(0.23,1.61} 


0.49 


0.11,2.77} 


253(205,307} 


177(121,230} 


8904448 


0.865981 


548.1 


0.41 


0.25,0.49} 


0.018(0.014,0.025} 


0.11(0.08,0.26} 


0.15 


0.09,0.76} 


171(158,192} 


11(6,15} 


8938628 


6.86219 


388.1 


0.22 


0.17,0.34} 


0.015(0.003,0.171} 


0.37(0.05,1.50} 


1.45 


0.15,4.05} 


127( 75,287} 


318(256,481} 


9451096 


1.25039 


106.7 


0.39 


0.25,0.65} 


0.069(0.019,0.283} 


0.49(0.13,1.33} 


0.61 


0.12,3.14} 


90( 59,144} 


66(42,107} 


9714358 


6.47418 


103.7 


0.27 


0.21,0.35} 


0.028(0.004,0.142} 


0.39(0.06,1.50} 


1.04 


0.15,3.91} 


65( 35,112} 


1252(1041,1558} 


9722737 


0.418528 


443.9 


0.55 


0.36,0.64} 


0.068(0.063,0.073} 


0.22(0.16,0.52} 


0.18 


0.09,0.92} 


230(225,236} 


2.4(1.6,2.8} 


9912977 


0.943916 


753.7 


0.23 


0.14,0.27} 


0.002(0.002,0.003} 


0.04(0.03,0.11} 


0.14 


0.08,0.66} 


105( 94,117} 


3.2(1.9,4.0} 


10095512 


6.01720 


472.6 


0.50 


0.37,0.71} 


0.185(0.072,0.579} 


0.88(0.22,2.15} 


0.78 


0.13,3.18} 


337(247,493} 


414(304,572} 


10226388 


0.660658 


934.9 


0.60 


0.39,0.72} 


0.124(0.101,0.150} 


0.35(0.23,0.83} 


0.24 


0.09,1.30} 


465(434,493} 


3.3(2.2,4.1} 


10319590 


21.3216 


247.1 


0.22 


0.10,0.62} 


0.013(0.001,0.642} 


0.34(0.04,2.05} 


1.08 


0.17,3.65} 


90( 34,329} 


4193(2175,9999} 


10613718 


1.175880 


88.1 


0.47 


0.30,0.72} 


0.136(0.063,0.449} 


0.75(0.26,1.73} 


0.75 


0.14,3.33} 


99( 76,147} 


80(52,121} 


10991989 


0.974478 


554.2 


0.54 


0.34,0.63} 


0.059(0.049,0.072} 


0.21(0.15,0.49} 


0.18 


0.09,0.92} 


256(239,274} 


11(7,13} 


11042923 


0.390164 


839.0 


0.40 


0.21,0.47} 


0.017(0.015,0.019} 


0.10(0.08,0.37} 


0.15 


0.09,1.37} 


223(213,230} 


< 1 


11968490 


1.078899 


253.2 


0.63 


0.43,0.80} 


0.333(0.287,0.387} 


0.88(0.55,1.69} 


0.52 


0.14,2.20} 


271(256,283} 


38(26,49} 



Note. — (1) The binary period is referenced to an epoch of BJD = 2454900; (2) Denned as sin 3 i (rip /(M 3 +M bin ) 2 , (3) See eq. J7J for the definition, (4) See eq. |lo| 
for the definition. The values in curly brackets represent the 10% lower- and 90% upper-limits on the probability distribution. The parameter values and uncertainties 
reported in this table are based on 10 s parameter draws for a 5-day binary, and scaled proportionally to /^,i n . 
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Table 3 

Candidate Triple-Star Systems Found in Kepler Data Base 



Source 


eccentricity 1 




r (3) 




i,,, 4 


v„, 5 


rms 6 


'Hongterm 


Rets. 








(degrees) 


(days) 


(d 


egrees) 


(degrees) 


(sec) 


years 




3228863 


U.UO 


U.UO,U. 1Z ) 


209(192,224} 


94( 63,123) 


45.4 


18.4,71.5) 


92( 13,139} 


51 


1600 




4647652 


U.J J 


U. 1U,U.44 J 


184( 42,340) 


459(113,644} 


44.9 


19.5,70.4} 


90( 21,160) 


35 


1400 




4909707 


U. J4 


U.J 1 ,U.OO J 


269( 17,359) 


420 ( 23,489) 


43.9 


24.2,63.7} 


88( 22,158) 


126 


305 




4940201 


U. 1 o 


nun o«; \ 

U. 1 1,U.ZJ ) 


163( 42,326) 


319(289,340) 


16.3 


9.2,21.4} 


54( 18,150} 


167 


41 




5039441 


42 


n i fi ft. l 

U. 1 o,U. J4 J 


187] 36,345) 


336( 48,619) 


45.4 


24.2,66.4) 


87( 21,159} 


39 


566 




5128972 


U. J J 


25 41 1 


101( 84,116) 


26( 7, 46) 


45.0 


18.4,71.4) 


86( 16,157) 


39 


1086 




5264818 


U.J / 


U. 1 J,U.J J j 


173{ 34,332) 


120( 23,270} 


41.4 


23.2,59.0) 


84( 22,154} 


62 


127 




5310387 


n 

U.JJ 


n ia n ai \ 

U.J4,U.Ol j 


161 ( 16,345) 


126( 12,194} 


45.7 


22.6,68.0} 


86( 21,161} 


20 


285 




5376552 


40 


35 45 ] 


167(161,175) 


302(296,309} 


44.3 


20.5,68.4) 


77( 16,171) 


39 


604 




5384802 


U.JO 


U.ZJ,U.40 ) 


171{ 30,334) 


103( 98,112} 


17.1 


9.3,23.4) 


84( 30,159} 


105 


29 


8 


5771589 


U.JU 


U.ZS,U. JJ J 


214( 38,329) 


75 ( 74, 76) 


31.4 


30.7,32.1} 


169(165,172} 


260 


3.2 


9 


6370665 


U.ZZ 


U.U /,U. JJ J 


92( 15,353) 


77{ 9,277) 


46.3 


23.1,67.7} 


68 ( 20,140} 


62 


240 




6525196 


U. JU 


U. ZO, U.JJ J 


285(233,310} 


187(127,200} 


28.0 


22.6,33.9) 


129( 84,147} 


29 


138 




6531485 


n aa 


U.JJ, U.OJ ) 


315(204,347) 


35( 33,35) 


37.8 


14.1,48.8) 


23( 8,175} 


68 


9.5 




6545018 


n 9^ 
u.zo 


n 1 £ n ia 1 

U. 10, U.JO ) 


150( 41,319) 


69( 67,71) 


21.8 


16.8,27.7) 


46( 23,63) 


109 


9 




7289157 


U.JO 


27 47 1 


161 ( 42,320) 


44(34,51) 


22.6 


15.3,29.7) 


68{ 9,172} 


73 


31 


9, 10 


7668648 


U. JO 


U.ZS,U.4Z J 


185( 40,327) 


29 ( 20, 36) 


36.8 


30.5,40.8) 


70( 59, 81} 


1193 


4 


9, 10 


7690843 


n 9 s 

U.Zj 


U.US,U.4Z J 


258( 48,334) 


44( 25,59) 


29.1 


17.1,42.2} 


101 ( 35,149} 


36 


19 




7837302 


n 1 ^ 

U. 10 


U.Uo,U.Zj ) 


247(175,319} 


353(302,397) 


14.7 


11.4,18.8) 


140( 15,169} 


120 


106 




7955301 




U.4J,U.4o j 


161 ( 36,326} 


187(186,188} 


31.6 


30.8,32.4) 


157(153,161} 


326 


8 


9 


8023317 


n 9^ 

U.Z J 


n ion 10 1 
u. is,u.zy j 


207 { 63,336} 


118( 92,145} 


53.0 


45.8,62.4) 


68 ( 52, 85} 


19 


62 




8043961 


n 9 s 

U.Z J 


U. 14, U.JJ ) 


192(167,212) 


398(363,425} 


34.6 


16.4,54.7} 


102( 11.172} 


50 


400 




8192840 


U.Oj 


U.JZ,U. 


173(160,185) 


569(544,595) 


45.0 


18.7,71.3} 


79( 24,164} 


59 


4108 




8386865 


U.JO 


f) T7 ft /Ifi l 
U.Z /,U.45 ) 


137(105,159} 


128(111,147) 


53.2 


33.1,74.0} 


120( 70,158) 


1 15 


187 




8394040 


U.Ol 


U.jU,U.O / ) 


123(113,131} 


296(288,305) 


43.8 


17.8,70.8) 


73( 19,159} 


96 


1088 




8719897 


0.24 


0.13,0.31 } 


291(267,317) 


90( 68,103) 


17.4 


( 9.2,25.2} 


98{ 29,151) 


51 


96 




8904448 


0.59 


0.50,0.66} 


135(125,143} 


443(431,454} 


40.1 


18.3,63.9) 


68( 12,166} 


32 


950 




8938628 


0.31 


0.26,0.35} 


282(221,327} 


339(314,348} 


17.4 


12.4,21.1 } 


133( 27.160} 


21 


60 




9451096 


0.24 


0.10,0.36} 


183( 53,313) 


60( 8,97) 


23.4 


11.9,37.1} 


91{ 33,150} 


19 


25 




9714358 


0.26 


0.20,0.32} 


154( 29,329) 


77 ( 76, 78) 


16.8 


13.8,20.8} 


134(120,149} 


131 


4.6 


9 


9722737 


0.22 


0.16,0.27} 


29 { 14,46} 


423 ( 5,441} 


45.1 


18.4,71.8} 


76( 14,161) 


48 


1290 




9912977 


0.31 


0.16,0.39} 


251(213,301} 


260(187,359} 


45.0 


18.7,71.2} 


103( 24,159) 


22 


1650 




10095512 


0.18 


0.12,0.23} 


67 ( 37,101) 


430( 19,457) 


13.6 


I 6.9,18.7) 


89{ 28,150} 


23 


100 




10226388 


0.32 


0.24,0.39} 


281(263,300} 


755(713,797) 


44.9 


18.4,71.9} 


80( 21,158} 


101 


3588 




10319590 


0.14 


0.05,0.32} 


182( 39,327} 


95( 82,111) 


10.4 


6.6,21.3} 


102( 11,171) 


470 


7.8 


9 


10613718 


0.18 


0.05,0.29} 


240(138,291} 


20( 7, 76} 


18.1 


9.7,29.0) 


121( 26,157) 


66 


18 




10991989 


0.30 


0.21,0.37} 


189(178,202) 


24( 6,546} 


43.0 


18.9,68.5) 


128 ( 21,165} 


82 


861 




11042923 


0.17 


0.09,0.25} 


41( 11,355) 


679(587,747} 


45.3 


19.6,71.2} 


92{ 25,162) 


57 


4950 




11968490 


0.40 


0.31,0.46} 


117(107,127} 


216(209,224) 


32.6 


16.2,48.9} 


57 { 29,128) 


43 


162 





Note. — (1) Orbital eccentricity of the triple system; (2) longitude of periastron (of the binary component) of the triple; (3) time of periastron 
passage of the triple; (4) mutual inclination angle between the orbital planes of the binary and triple; (5) ang le between the triplets periapse and 
the plane of the binary (see Fig.[6j - v m runs between 0° and 180° because of the way it appears in eq. j9J; (6) rms scatter of the O — C p oint s 
about the best- fitting model; (7) timescale for longer-term perturbations in the triple system calculated here simply as P^ p /Pb\ n (see eq. 1 131; 
(8) Fabrycky (2010); (9) Slawson et al. (2011); (10) Carter et al. (2013). The values in curly brackets represent the 10% lower- and 90% 
upper-limits on the probability distribution. The parameter values and uncertainties reported in this table are based on 10 8 parameter draws for 
a 5-day binary, and scaled proportionally to i^,i n . 



